#ifndef __FEYNRULE_H__
#define __FEYNRULES_H__
#include <ginac/ginac.h>
#include <bubblesii/pfrac.hpp> // header with prop class
#include "field.h"
#include "symconst.h"
using namespace cnst;
using namespace GiNaC;
/// Templates for expr generation and
/// reduction number of graph vertices
/// Base classes for feynman rules functions
class gauge_boson
{
 public:
  int t;
  ex momentum;
  int Fr;
  int To;
  field from_part;
  field to_part;
  ex expr;
  gauge_boson(int id1,int id2,ex mom,ex mass,field from_part_,field to_part_)
    { 
      t = 0;
      Fr = id1;
      To = id2;
      from_part = from_part_;
      to_part = to_part_;
      momentum = mom;
      expr = bubblesii::prop(mass,1,momentum)*lorentz_g(mu(id2).toggle_variance(),mu(id1).toggle_variance());
    }
};

// gauge boson propagator inarbitrary gauge
class gauge_boson_xi
{
 public:
  int t;
  ex momentum;
  int Fr;
  int To;
  field from_part;
  field to_part;
  ex expr;
  gauge_boson_xi(int id1,int id2,ex mom,field from_part_,field to_part_,ex mass,ex xi_)
    { 
      t = 0;
      Fr = id1;
      To = id2;
      from_part = from_part_;
      to_part = to_part_;
      momentum = mom;
      expr = -I*bubblesii::prop(mass,1,momentum)*( lorentz_g(mu(id2).toggle_variance(),mu(id1).toggle_variance())
                                                   -(1-xi_)*indexed(momentum,mu(id1).toggle_variance())*indexed(momentum,mu(id2).toggle_variance())*bubblesii::prop(xi_*mass,1,momentum));
    }
};

class scalar
{
 public:
  int t;
  ex momentum;
  int Fr;
  int To;
  field from_part;
  field to_part;
  ex expr;
  scalar(int id1,int id2,ex mom,ex mass,field from_part_,field to_part_)
    { 
      t = -1;
      Fr = id1;
      To = id2;
      from_part = from_part_;
      to_part = to_part_;
      momentum = mom;
      varidx tmp_idx8(symbol("dnu8"),4);
      //expr with indices
      ///!       expr = 1/(indexed(mom,tmp_idx8)*indexed(mom,tmp_idx8.toggle_variance()) - pow(mass,2));      
      expr = bubblesii::prop(mass,1,mom);
      //expr without indices
      //expr = 1/(mom*mom - pow(mass,2));      
    }
};

class fermion
{
 public:
  int t;
  ex momentum;
  int Fr;
  int To;
  field from_part;
  field to_part;
  ex expr;
  fermion(int id1,int id2,ex mom,ex mass,field from_part_,field to_part_)
    { 
      t = 1;
      Fr = id1;
      To = id2;
      from_part = from_part_;
      to_part = to_part_;
      momentum = mom;
      expr = (dirac_slash(mom,4-2*ep,id2+1) + mass*dirac_ONE())*bubblesii::prop(mass,1,mom);
    }
};


/// Propagators
template<field p1,field p2>
  class prop{
  };

template<>
class prop<e,E> : public fermion
{
 public:
 prop(int id1,int id2,ex mom) : fermion(id1,id2,mom,0/*me*/,E,e)
    {
      //varidx tmp_idx1(symbol("dnu1"),4);
      //ex e = (dirac_slash(mom,4) + sim*dirac_ONE())/(pow(mom,2) - pow(sim,2));
      //      expr = (dirac_slash(mom,4) + me*dirac_ONE())/(indexed(mom,tmp_idx1)*indexed(mom,tmp_idx1.toggle_variance()) - pow(me,2));
      //  ex e = (dirac_slash(mom,4) + sim*dirac_ONE())/((dirac_slash(mom,4) - sim*dirac_ONE())*(dirac_slash(mom,4) + sim*dirac_ONE()));
    }
};

template<>
class prop<A,A> : public gauge_boson
{
 public:
 prop(int id1,int id2,ex momentum) : gauge_boson(id1,id2,momentum,0,A,A)
    {
      //      varidx tmp_idx1(symbol("dnu1"),4);
      // expr = lorentz_g(mu(id2).toggle_variance(),mu(id1).toggle_variance())/
      // (indexed(momentum,tmp_idx1)*indexed(momentum,tmp_idx1.toggle_variance())-pow(m0,2));
    }
};

// propagator in arbitrary gauge for Machacek
template<>
class prop<A_xi,A_xi> : public gauge_boson_xi
{
 public:
 prop(int id1,int id2,ex momentum) : gauge_boson_xi(id1,id2,momentum,A,A,0,xi) //? mb A_xi?
    {
    }
};


/*
template<>
class prop<phi,phi> : public scalar
{
 public:
 prop(int id1,int id2,ex momentum) : scalar(id1,id2,momentum,mS,phi,phi)
    {
    }
    };
*/
      
                      // bubblesii prop variant
template<>
class prop<phi,phi> : public scalar
{
 public:
 prop(int id1,int id2,ex momentum) : scalar(id1,id2,momentum,0,phi,phi)
    {
    }
};

///Vertices 3 part   t=momentum power
template<field f1,field f2,field f3>
  class vert3{
};

template<>
class vert3< E,e,A >
{
 public:
  int t;
  ex expr;
   vert3(int id,ex mom1,ex mom2,ex mom3)  
    {
      expr = dirac_gamma(mu(id),id+1);
      t = 0;
    }
};

template<>
class vert3< phi,phi,A >
{
 public:
  int t;
  ex expr;
   vert3(int id,ex mom1,ex mom2,ex mom3)  
    {
      expr = -I*indexed(mom1-mom2,mu(id));
      t = 1;
    }
};

template<>
class vert3< phi,phi,A_xi >
{
 public:
  int t;
  ex expr;
   vert3(int id,ex mom1,ex mom2,ex mom3)  
    {
      expr = -I*indexed(mom1-mom2,mu(id));
      t = 1;
    }
};
  
///Vertices 4 part
template<field f1,field f2,field f3,field f4>
  class vert4{
};

template<>
class vert4< phi,phi,phi,phi >
{
 public:
  int t;
  ex expr;
  vert4(int id,ex mom1,ex mom2,ex mom3,ex mom4)  
    {
      expr = -I*lB;
      t = 0;
    }
};

#endif // __FEYNRULES_H__
